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ABSTRACT 

The inner two planets around the 55 Cancri were found to be trapped in the 3:1 mean 
motion resonance. In this paper, we study the dynamics of this extra-solar planetary 
system. Our numerical investigation confirms the existence of the 3:1 resonance and 
implies a complex orbital motion. Different stable motion types, with and without 
the apsidal corotation, are found. Due to the high eccentricities in this system, we 
apply a semi-analytical method based on a new expansion of the Hamiltonian of the 
planar three-body problem in the discussion. We analyse the occurrence of the apsidal 
corotation in this mean motion resonance and its influence on the stability of the 
system. 
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1 INTRODUCTION 

Over 100 extra-solar planets have been found (see e.g. 
http://www.obspm.fr/encycl/encycl.html by J.Schneider). 
The high eccentricities, close distances to the central stars 
and the heavy masses, show the wide variety of planetary 
system that are quite different from our solar system. 

Among these "exoplanets" , 28 are located around the 
main sequence stars in 13 "multiple planet systems", and 
some of them appear to be in Mean Motion Resonances 
(hereafter MMR). Recently, the two inner planets around 
the 55 Cancri were found to be in a possible 3:1 MMR (Ji et 
al. 2003). In numerical simulations, the resonant angles are 
found to librate around certain values. At the same time, 
the difference between the two periastrons is found to be 
locked to a fixed value, that is, the two planets precess at 
the same rate. 

It is commonly accepted that these extra-solar planets 
are not formed in situ but have suffered orbital migrations. 
Different migration speeds of planets may lead to a vari- 
ation of mutual distance between two neighboring planets 
and eventually a capture into the current resonant configu- 
rations (Kley 2000; Kley 2003; Nelson & Papaloizou 2002). 
The dynamics of these resonances may contain important in- 
formation of the capturing process and therefore give valu- 
able hints to the planet formation and the evolution of a 
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planetary disk. As a result, planetary systems with reso- 
nance require special attention. Moreover, the possibility of 
an Earth-type (habitable) planet in the 55 Cancri system 
has been discussed (Cuntz et al. 2003). The existence of a 
habitable planet in a system depends not only on the cur- 
rent dynamical features of the system but also on the his- 
tory it has experienced. Therefore it is important to study 
the now-known planets before an Earth-like planet was ob- 
served someday. All these points make this planetary system 
interesting. 

In this paper, we confirm the possibility that the two 
planets are in a 3:1 MMR by numerical integrations and 
discuss the possible configurations this system would take 
(section 2), then with a new expansion of the Hamiltonian 
for a planar three-body problem (section 3), we analyse the 
dynamics of this system. We will discuss the occurrence of 
the apsidal corotation and its influence on the stability of 
this system (section 4). Finally the conclusions and discus- 
sions are given in section 5. 



2 NUMERICAL SIMULATIONS 

2.1 Numerical model and initial set-up 

The 55 Cancri system is a system with three planets. We 
adopt the orbital solution given by Fischer et al. (2003), as 
being listed in Table 1. To numerically simulate this sys- 
tem, we use a symplectic integrator (Wisdom & Holman 
1991; Mikkola & Palmer 2000), which allows us to follow 
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Table 1. The orbital solution from Fischer et al (2003). The 
mass of the central star is 1.O3M0, and in this paper we adopt 
the planet masses by assuming sini = 1. 



Parameter 


Companion B 


Companion C 


Companion D 


Msini (Mj) 


0.83 


0.20 


3.69 


P (days) 


14.65 


44.27 


4780 


a (au) 


0.115 


0.241 


5.461 


e 


0.03 


0.41 


0.28 



the orbital evolution of each body and simultaneously the 
stability of this orbit indicated by the Lyapunov Character 
Indicator (LCI) (Froeschle 1984). Each planet orbit has an 
LCI and we choose the largest one among the three to in- 
dicate the stability of this planetary system in this paper. 
The time step is set to be 0.3 d, which is about 2% of the 
orbital period of the innermost planet. 

We adopt the masses, semimajor axes and eccentricities 
of planets the values listed in Table 1. The initial orbital 
inclination of companion D is set to be zero, while the com- 
panion B and C have initial inclinations of 10 degrees. 
Other angles (the ascending node, the periastron and the 
mean anomaly) are randomly generated from [0, 2ir). Four 
hundred simulations, with different initial conditions, are in- 
tegrated up to 10 6 years. 

Ignoring the companion D, we first integrate a system 
composed of the central star and the companion B and C. 
Then, retaining the initial conditions of these three bodies 
and adding the companion D on an initial orbit with the 
given a, e, i and randomly selected angle elements, we inte- 
grate the four-body system. As the simulations show, the 
companion D has nearly no influence on the motion of the 
inner planets, although it is much heavier. This is due to 
its large semimajor axis (5.461 au). However, in order to get 
reliable simulations for the real system, we report here the 
numerical results from the four-body model. 

During the integrations, if the distance between any two 
of the planets became smaller than half of the criterion for 
the "Hill stability" (Gladman 1993), the system was consid- 
ered become collapsed and the simulation was terminated. 

Generally, the motion of a planet is dominated by the 
central star and its orbit is a conic section with small de- 
viations due to the gravitational perturbation from other 
planet(s). This perturbation can be described by the disturb- 
ing function, which can be expanded in terms of the orbital 
elements. We use a, e, i, w, fl, A to denote the semi-major 
axis, eccentricity, inclination, longitude of periastron, longi- 
tude of ascending node, and mean longitude, respectively. 
The disturbing function for a planet with mass mi and or- 
bital elements (ai,ei,ii,ti7i,fli,Ai), perturbed by another 
planet (indicated by subscripts '2'), can be written as 

R = (12 ^ F(a 1 ,a 2 , ei, e 2 , ii, 12) cos ip. 

Here \i2 = Qm2 and <p is a linear combination of angles 

V = JlAl + 32X2 + ja^l + + J5^1 + J6^2 

where ji , 32 , • • • , ja are integers satisfying ji +32 H Vje = . 

Particularly, when the orbital periods satisfy P2/P1 ~ 
3, that is, Ai — 3A2 ~ (dots denote the time derivation), 
we have <P(j 1 =i,j 2 =-3) ~ since generally «7i, W2, £2i, &2 *C 



Figure 2. A typical temporal evolution of the semimajor axes 
ai , C12 in a stable system. 

Ai, A2. Then in the disturbing function 7?, those terms con- 
taining Ai — 3A2 can not be eliminated by the averaging 
technique and become the leading terms. And if the in- 
clinations 11,12 are not too large, the leading terms in R 
are of 0(e\), 0(e%), 0(e\e2) and the corresponding angle ip, 
now called resonant arguments, are 81 = Ai — 3A2 + 2zvi, 
62 = Ai — 3A2 + 2ci72 ,83 = Ai — 3A2 + vj\ + -UJ2 ■ Therefore the 
libration of any one of #1,2,3 indicates the two planets are in 
an {eccentricity-type) 3:1 MMR. We note that no more than 
two of them are linear independent and 83 = (8\ +0 2 )/2 [see 
for example Murray & Dermott (1999) for more details]. 

Beside these resonant arguments, the relative apsidal 
longitude Avj — w\ — zu 2 is also a critical argument often 
to be discussed. We call the libration of Azu in an MMR 
the apsidal corotation. In a 3:1 MMR, if both 81 and 82 
librate with small amplitudes, Avj = \(8\ — 82) will also 
necessarily librate. In this sense, the apsidal corotation is not 
independent and should be different from the term of apsidal 
resonance, which is in the context of a secular perturbation. 

Hereafter we use subscripts 1, 2 to label the orbital ele- 
ments of companion B and C respectively. 

2.2 Numerical results 

About one third (133 out of 400) of the simulations collapse 
during the integrations, and all the remainders survive for 
the 10 6 yr integration. Most of the survivors have very short 
(< 10 2 yr) e-folding time T e (= 1/LCI, the reciprocal of 
LCI), while 38 of them have T e ^ 10 3 yr, and they are m 
this sense regarded as stable. In these stable systems, all the 
three planets have final inclinations ^ 0°.l, but those un- 
stable (T e < 10 yr) systems may have inclinations of com- 
panion B and C as high as ~ 20°, that is, the stable system 
prefers to retain to be coplanar. All the stable systems are 
found to be associated with the 3:1 MMR between the two 
inner planets. 

According to different configurations of #1,2,3 and Aw, 
the stable systems can be divided into three groups, with 
representative cases illustrated in Fig. 1. In all the cases, 
the companion D is in a steady motion (we don't illustrate 
this here) and, the semi-major axes of companion B and C 
vibrate with very small variations (an example is shown in 
Fig. 2). 

Case a (Fig. 1) shows a strong resonance with 81, 82 
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Figure 1. The temporal evolution of the critical angles 61, 02, 03, Azu, and the eccentricities of planets ei,e2 (in the bottom panel, the 
dashed curves indicate e^)- Case a, b, and c are from different initial conditions. 



and 83 librating around 215°, 75° and 325°, respectively. 
Simultaneously, the Azu also librates around 250° . The well 
known symmetric configurations of Azu librating around 0° 
(alignment) or 180° (anti-alignment) imply the conjunctions 
of the two planets always occur at a certain position and this 
attributes to the stability of the system (Lee & Peale 2002) . 
However the asymmetric configuration with Azu librating 
around angles other than 0° or 180°, has been predicted by 
Beauge, Ferraz-Mello & Michtchenko (2003). 

Case b differs from case a in the behavior of 8\ , 83 and 
Azu, which circulate now. The 82 still librates around 80° 
with an amplitude of 40°. In fact, when 82 librates with a 
small amplitude, we may assume 62 ~ const., that is, 82 ~ 0. 
From 82 — Ai — 3A 2 + 2rij 2 ~ we can derive 

ri2 — ZU2 ^ 1 

m — Z02 3 ' 

where ni j2 = A ii2 are the mean motions of planets. In a 
reference frame corotating with the periastron of the second 
planet, the relative mean motions are n' k = n k ~ ZU2 (k = 
1,2), and from n'2/n'i w 1/3 we know that in this reference 
frame the conjunctions of planets always happen at a certain 
direction. Meanwhile, given n k >• ZU2, we have 112/ni « 
n'2/n'i and there is still a 3:1 relation between the two orbital 
periods. In the 55 Cancri system, our numerical results show 
ni,2 are larger than zb2 by at least two order of magnitudes. 

Although the numerical simulations give more examples 
with configuration as case a than case b (28 vs 7, see Table 
2), we notice that among the most stable sets with e-folding 
time T e ~ 10 6 yr (the same order as the integration time), 
there are two sets having configurations as case b but only 
one as case a. This prevents us from concluding which type, 
case a or case b, is more stable. 

In case c, 81,83 and Azu circulate, while 82 librates 
around 180° with a quite large amplitude of 160°. As a re- 
sult, it is less stable than the above two cases. The only three 



Table 2. The distribution of T e of the stable systems. The num- 
bers with subscripts a, b or c are the sums of systems with dif- 
ferent configurations as shown in Fig. 1. 



Te (yr) 


~ 10 3 


~ 10 4 ~ 10 5 ~ 10 6 


Number 


12 a + 2 b + 3 C 


8a + lb 7 a + 2 b l a + 2 b 



examples with this configuration in our simulations all have 
e-folding times < 5 x 10 3 yr. 

Finally we summarize the e-folding times and the mo- 
tion configurations of the 38 stable systems in Table 2. 

Besides the angles #1,2,3 and zu, we notice that the or- 
bital eccentricities ei,e 2 also have remarkably different be- 
haviors in these three cases. Another noteworthy feature of 
numerical results is that the system has a symmetry over 
angles 8i and Azu, that is, when is a libration center of 
8i or Azu, the argument (360° — 0) is too. For example, we 
observe 82 librates around ~ 80° and ~ 280° with nearly 
the same ratio in our simulations. 



2.3 Some extended calculations 

To integrate all the surviving systems to 10 7 yr or even 
longer is out of our current computational capability. How- 
ever, to explore further the stabilities of different configura- 
tions, we integrate the three examples shown in Fig. 1 further 
to 3 x 10 7 years. Both the case a and b last this time span 
and keep the configurations all the time. But case c loses 
the stability when the inclinations of companion B and C 
increase to ~ 5° after ~lx 10 7 years. As in those unstable 
examples, this inclination increasing is due to a sequence of 
close encounters between the two planets. On the contrary, 
the mechanism of the 3:1 MMR can protect those systems 
as case a and b from catastrophic close encounters. 

Other 5 examples randomly selected from the 17 sys- 
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Table 3. The percentages of the survival and stable systems in 
the simulations with different initial inclinations (i® = 12)- 



»1,2 ( dc §-) 


10~ 5 


2.5 


5 


10 


15 


20 


25 


p 

r survival 


66 


62 


52 


35 


23 


23 


21 


^stable 


10 


10 


10 


10 


6 


5 


6 



tems with T e ~ 10 3 yr are integrated to 1 x 10 7 yr. Again, ex- 
cept the one having motion configuration as case c, they all 
survive and have new e-folding time T e > 10 3 yr. [The LCI 
will converge to a certain value as the integrating time ap- 
proaches to infinity, so T e of an orbit would change when it's 
integrated to a longer time span, see for example (Froeschle 
1984).] As for those systems with larger T e , they are more 
stable and should have longer surviving lives. 

So far, our discussions are all based on the small initial 
inclinations of the two inner planets. We find that in stable 
systems the inner two planets are locked to the 3:1 MMR and 
their small orbital inclinations can be kept to at least 10 6 yr. 
While in those unstable systems, chaotic close encounters of 
planets could cause inclinations to increase. Furthermore, we 
also numerically simulated hundreds of examples with larger 
initial inclinations. In these systems, except the inclinations 
of the two inner planets are set to higher values, all the 
orbital elements are initialized in the same way as before, 
and we integrate these systems to 10 6 yr. A brief summary 
is given in Table 3. 

We see the survive probability drops from 66% to 21% 
as the initial inclinations increase from 10" 5 to 25 degrees, 
but at the same time the percentage of stable systems even 
does not decrease until i? j2 = 10° before it attains a value 
of ~ 5%. This implies that the stability of the system is 
not so sensitive to the inclinations, and if smaller than 10 
degrees, they have nearly no influence on the stability of the 
system. On the other hand, that all the stable examples are 
associated with the 3:1 MMR reveals again the importance 
of this resonance in stabilizing this planetary system, even 
when the orbits are far from coplanar. We don't give more 
analysis to the problem of higher inclinations here, but it 
deserves a detailed discussion in future . 



3 ANALYTICAL MODEL 

To understand the dynamics of the 3:1 MMR in the 55 Can- 
cri, it is useful to describe it with an analytical model. As 
we mentioned above, the companion D has a very weak in- 
fluence on the motion of the inner planets. In addition, the 
numerical simulations show it is mainly the 3:1 MMR that 
is responsible for the stability of the system, in spite of the 
higher orbital inclinations. So, a planar three-body prob- 
lem is a reasonable model for describing this planetary sys- 
tem, especially when we mainly focus on the 3:1 MMR. Of 
course, considering the possible substantial difference be- 
tween the coplanar and non-coplanar motion, a more realis- 
tic and reliable model is still needed. On the other hand, the 
usual analytical perturbation models, such as the Laplace 
expansions, are not convergent for high eccentricities, and 
expansions with low-order truncations may yield imprecise 
results (Ferraz-Mello 1994). Recently, a new expansion of the 



Hamiltonian suitable for high-eccentricities has been devel- 
oped by Beauge and Michtchenko (2003). We will apply this 
expansion to the system the 55 Cancri. 

3.1 The Hamiltonian 

Suppose three bodies with mass Mo, mi, mi orbiting around 
their common center of mass. Mo 3> mi,2 is the central star 
and mi, 2 > are the planets. Near the 3 : 1 MMR, a set of 
canonical variables are defined as: 



Ai; 




Ji 


= Li + i(/i+/ 2 ) 


A2; 




J2 


= La-|(/i+/2) 


3, 


I- 


— vji ; I\ 


= Li-Gi (1) 


C"2 = — 


¥• 


— TO2; I2 




where Li = 


m' iy 


[Juai and G; = 


= Li yjl — ef are Delaunay 



variables, Oi,e», A;, to; are the orbital elements of the ith 
planet, [n — Q(rm+Mo), Q is the gravitational constant, and 
/ _ m.Mo j s -j-jjg re duced mass. These canonical variables 
define a four degree-of- freedom system. And let's bear in 
mind 9i = — 2ct; (i — 1, 2), #3 = — ai — 02 and Azu = cr 2 — ai. 

After considering the conservation of angular momen- 
tum and averaging over the fast-angle, we get two integrals 
of motion J aU m = Ji + J2, JdU = Ji — Ji- At the same time, 
the problem is reduced to a two degree-of-freedom system in 
canonical variables (i;,cr;). The resulting Hamiltonian reads 

^=-e4#+^ ( 2 ) 

1=1 1 

where Hi is the truncated disturbing function, 

3 j max fe max U max i max 

Hl = ""a™ 2 EEEE 12 i?nj, fc ,.,i(a-a ) n 

n=0 j=0 fc = 11=0 i = -i m ax 

xe{e2 cos(2ti(Ji + l(a 2 — fi))- (3) 

Here a = ai /a 2 is the ratio between the two semimajor 
axes and a is its value in the exact MMR. The coeffi- 
cients Rn,j,k,u,l can be determined beforehand (Beauge & 
Michtchenko 2003). Considering the high value of eccentric- 
ities, we adopt the upper limits of sums tt m ax = (max = 12 
and jmax = fcmax = 15 in this paper. For convenience, the 
Hamiltonian has been expressed explicitly in (oj, e;) instead 
of (Ii, Ji), therefore the integrals J sum and Jdif don't appear 
explicitly but exist as constrains. We note the total energy 
H is an integral too. Furthermore, it is easy to prove the 
symmetry buried in the Hamiltonian: 

#(<ti, <r 2 ) = H(-ai,-a 2 ) = ff(180° + ai, 180° + a 2 ). 

3.2 The reliability of the Hamiltonian 

We can calculate planet orbits by applying this Hamiltonian. 
Since it is an averaged system, the time step can now be 
much longer than that used in our symplectic integrator. 
A comparison between this Hamiltonian expansion and the 
direct numerical integration is shown in Fig. 3, in which we 
choose the most unstable case of Fig. lc as an example. We 
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Figure 3. The temporal evolution of 9 2 , Aro and eccentricities 
ei,2 from the same initial conditions as the ones of Fig. lc. Squares 
show the results of a numerical integration of the exact dynami- 
cal equations for a 4-body model, while the curves are from the 
Hamiltonian equations for a 3-body model. 

get the same initial conditions and recalculate the orbits 
from the Hamiltonian equations with an average (variable) 
time step of 60 d, which is four times of the inner planet's 
orbital period. 

The agreements in Fig. 3, not only imply the reliability 
and suitability of this Hamiltonian, but also reflect that the 
third planet can be ignored when studying the dynamics of 
the inner two planets. 



4 THE APSIDAL COROTATION AND THE 
STABILITY OF THE SYSTEM 

The Hamiltonian introduced above has two degrees of free- 
dom. Although it is still non-integrable, it can help us to 
understand the dynamical properties of the system. In this 
section, we will use it to find out why and when the system 
follows the different types of resonant evolutions as shown in 
Fig. 1, and then via the surface of section technique we study 
the stability of the systems with different configurations. 

4.1 The eccentricities and the resonant angles 

Since J sum and J di f are integrals, Ji = |(J sum + Jdif) 
and J 2 = i (J sum — Jdif) are invariables too. They con- 
strain the variation of the planetary semi-major axes by 
3Ji + J2 = 3m' lx //iiai + m^^/Jl^a^ = const. However, when 
we discuss planets trapped in an MMR, we may simply as- 
sume the semi-major axes are constants, otherwise the reso- 
nance would be destroyed. In fact, the numerical simulations 



Figure 4. The possible variation of eccentricities of the planets. 
The dotted curve corresponds to J sum = 3.951733 X 10~ 4 . Results 
from the numerical simulations of Fig. la (plus), b (black dot) and 
c (open circle) are plotted to compare. 

in section 2 also support this assumption (Fig. 2). Similarly, 
simple calculations show that 

Jsum = mi \J //iOi(l - ef) + m' 2 \J u 2 a 2 (l - e%). (4) 

This confines the variation of the eccentricities when we 
fix the semi-major axes a\ and a 2 . Besides this, the vari- 
ation range of the eccentricities is also bounded by other 
constrains, e.g., the total energy of the system. 

We assume the two semi-major axes are fixed, 01 = 
0.115 au and a 2 = 0.241 au. Then Eq. (4) defines a family 
of curves on the (ei, e2)-plane (Fig. 4). Each curve in Fig. 4 
corresponds to a definite value of the integral J SU m. The 
orbital elements in Table 1 give J SU m = 3.591733 x 10~ 4 , so 
that the eccentricities e\ , e 2 of orbits beginning from these 
initial conditions are restricted on the corresponding curve 
during the dynamical evolution. 

After assuming 01,2 fixed, for certain ei,2 values, we 
can calculate the energy level curves of the Hamiltonian on 
the (<7i, <j 2 )-pl&ne. Fig. 5 shows such contours. At the first 
glance, we note the symmetry over both <j\ and a 2 . This 
has been observed in our numerical simulations (recalling 
0i = -2oi). 

Since the total energy H is an integral, a system start- 
ing from a point on any curve in Fig. 5 can never jump to 
other curve with a different energy. On the other hand, the 
variation of eccentricities (Fig. 4) is constrained by the in- 
tegral J S um through Eq. (4) . Therefore, as the eccentricities 
vary following Eq. (4), an energy curve on the (d, o 2 )-p\anc 
will extend to an energy surface (e.g. a cylinder) and the 
system will evolve on this surface. To illustrate this, we fix 
di = 0.115 au and a 2 — 0.241 au, set ei = 0.03 (and se- 
quently 0.10, 0.15, 0.20), evaluate the corresponding e 2 from 
Eq. (4) with J sum = 3.591733 x 10~ 4 , then calculate contours 
for different pairs of (ei,e 2 ). Finally, we get a series of sec- 
tion of the energy surface in Fig. 5. In this way, we can follow 
the dynamical evolution of a system starting from definite 
initial conditions. 

For example, the thick curves in Fig. 5, representing an 
energy H — —3.957679 x 10~ 3 , keep the closed oblate shape 
when ei evolves from the initial value 0.03 to 0.20. The 
oblate profile of the curve implies that the angle ai (9i) has 
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Figure 6. The contour of Hamiltonian on the [ei, Ara)-plane. 
Some Hamiltonian values are labelled. Dots and squares represent 
the case a and b in Fig. 1 from numerical integrations. 



Figure 5. The energy level curves of the Hamiltonian on the 
(ci, <T2)-plane. A gray level code is used to measure energy val- 
ues, the darker regions indicate smaller energy value, and vice 
versa, a, b, c, and d are the case for e\ = 0.03, 0.10, 0.15 and 0.20 
respectively, with obtained from Eq. (4). 

a larger reachable range than a 2 (02)during the temporal 
evolution of a system with initial conditions located on this 
curve. The angle o\ could vary in a range of ~ (160°, 400°), 
while (72 is bounded between 20° and 60°. Correspondingly, 
the resonant angle Q\ and #2 could vary in ranges with widths 
of 480° (circulation) and 80° (libration), respectively. This 
is just the situation of case b in Fig. 1. Another example 
indicated by a dashed curve with H = -3.957722 x 10~ 3 
exhibits the possibility of circulating of both 0i and 02. 

We denote the possible variation of Oi in Fig. 5 by Acrj. 
According to the foregoing analysis, A 02 < Ao\. If Aoi <^ 
180°, both 0i and 8 2 would librate (A0 2 < A0i <, 360°) 
and the system could run as Fig. la. It could run as Fig. lb, 
if Acxi <; 180° while 02 is still bounded in a narrow range, 
say Acr 2 <, 90°. And if 90° <, Aa 2 <, 180°, Fig. lc could be 
followed. So the measures of such areas in the whole (<ti, 02)- 
plane tell the probabilities of different types of motion. Cal- 
culations show the proportions of these regions in the whole 
plane, i.e. the probabilities of motion type a, b and c, ap- 
proximatively are 10%, 16% and 4.5%, respectively. They 
are coarse estimates because we don't count the stability of 
the a possible motion. After considering the stability, these 
estimates must decrease in different extents. 

In a word, given ai,ei, the dynamical evolution of a 
system is determined by the initial values of Oi (i = 1, 2), 
which consequently correspond to different energy H. 

4.2 The occurrence of apsidal corotation 

When two planets are in a 3:1 MMR, at least one of the three 
resonant angles librates around a definite value. Since 02 
generally has a smaller libration amplitude than ui, we can 
further simplify the Hamiltonian by safely setting 02 = 0%. 
To determine a®, we can calculate the stable stationary so- 
lutions of the Hamiltonian equations (Beauge et al 2003), 
or we can get an estimate from the position of the Hamil- 
tonian extremum in Fig. 5. But in this paper, we simply 



adopt 02 = 320° (02 = 80°) from the numerical results. Re- 
calling the constant semimajor axes and the constrain on 
eccentricities, the Hamiltonian can be explicitly expressed 
as H — H(a2, a}{, e 2 , ei(J sum , £2), 02, Acer), with only two 
"free" variables e 2 and Aw. We calculate the contour of such 
Hamiltonian on the [e 2 , Am)-plane and show it in Fig. 6. 
The contour of the velocity of Aw is also calculated and 
shown in Fig. 7. 

A system approximately satisfying the assumptions 
mentioned above will evolve along one of the curves in Fig. 6. 
For the inner two planets in the 55 Cancri, the initial orbital 
elements listed in Table 1 (e? = 0.03, e° = 0.41) correspond 
to the right boundary of the box in Fig. 6. Two different fates 
of a system starting from these initial conditions are distin- 
guished: one is characterized by a libration of Aw around 
~ 250° and the other a circulation of Aw. In the latter case, 
the eccentricity e2 can reach a smaller value (correspond- 
ingly ei a higher value). These agree pretty well with the 
directly numerical integrations. Such agreements are shown 
in Fig. 6 by a comparing plots of the numerical results. 

Fig. 7 shows the velocity of Aw. Through this figure, 
we can explain how Aw evolves. Along the right bound- 
ary of the box, the absolute value of velocity could be 
very large, so that Aw increases or decreases quickly. Let 
us write the perturbing part of the Hamiltonian as Hi — 
X)j™q X Sfe™jT ^i. fee i e 2i where Rj,u contains the coefficients 
Rn,j,k,u,i and sums over n, u, I in Eq.(3). By careful calcula- 
tions, we obtain 

j=l k=0 

-^EE M ^t (5) 

2 j=0 k=l 

Obviously, when e; <g 1, the dominating term on the right 
hand side has the order ~ — . Near the right boundary of the 
box in Fig. 7, we have ei ~ 0.03 1. That's the reason why 
Aw has a quick change when ei approaches its minimum 
value as shown in Fig. lb,c. 

The numerical simulations reveal that the quick varying 
of Aw in Fig. lb,c is caused by a quick varying of w\ when 
ei — » 0. In fact, when an orbit is nearly circular (e ~ 0), 
it would be relatively easy to change its direction (w), thus 
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Figure 7. The contour of dAzu/dt. Some velocity values are la- 
belled and the scale is 360° yr _1 . Dashed curves indicate the zero 
velocity curve. Dots and squares are the same as in Fig. 6. 



Figure 8. The initial conditions of the 38 stable systems. Open 
squares, plus and solid squares indicate the systems evolving as 
case a, b and c in Fig. 1. The dotted lines give the positions of 
Azu = 160°, 330° and the solid line is a° = Azu + a° = 320°. 



the two planets could loose the locking of periastrons if they 
were not in a strong coupling. On the other hand, whether 
they can strongly couple with each other depends also on 
how close they can approach each other, that is, how large 
their eccentricities (especially e 2 ) are. Bearing this in mind 
we can understand why in case a e 2 has a higher lower-limit 
than in cases b and c. Of course, a small e 2 can also cause 
a quick variation of Azu, but this situation does not happen 
because the energy integral prevents e 2 from approaching 
very small value. 

From Fig. 6 and Fig. 7, we see different behaviors of 
Azu (circulation or libration) have different varying direc- 
tions along the right edge of the box. By numerically solv- 
ing the equation dAzu/dt = 0, we find two zero velocity 
points on the boundary. They define the boundaries of the 
initial conditions leading to a librating or circulating of Azu. 
The calculations show that the apsidal corotation happens 
if Azu e (f60°,330°), while if Azu G (-30°, 160°), Azu 
circulates. 

The case c has a Hamiltonian value a little further from 
the extremum in Fig. 5, so that <r 2 (f? 2 ) has a large amplitude 
of libration and can not be assumed as a constant. As a 
result, we can not find a corresponding curve for it in Fig. 6. 
However, if we draw a series of Hamiltonian contour with 
different value of cr 2 (adding another dimension to Fig. 6), 
we can also understand its dynamical evolution. 

Fig. 8 summarizes the initial conditions of the 38 stable 
examples in the numerical simulations. For clarity, we have 
mapped the symmetric configurations to one regime accord- 
ing to the symmetry of the Hamiltonian. Except the three 
examples as case c, all the points for stable systems gather 
along the line of a° = 320°, and points for librating and 
circulating Azu (case a and b) are separated by the dotted 
lines of Azu = 160° , 330° . These agree very well with the 
above analysis with the Hamiltonian. 

We have given an estimate of the probability of follow- 
ing different ways as case a, b and c. The analyses here tell 
that all the initial conditions leading to stable configurations 
gather in a narrow strip along the line in Fig. 8. That is, for 
systems having given ai_ 2 , ei, 2 as in Table 1, only those with 
initial conditions satisfying a° « 320° and a° + Azu w 320° 
can survive for long time. 



4.3 Surfaces of section 

We know now the initial conditions determine whether the 
apsidal corotation happens or not, but whether a system 
would follow a definite evolution type also depends on the 
stability of an evolution type. The stability of an orbit 
in the phase space can be revealed by its projection on 
the surface of section. We present in Fig. 9 the sections of 
(7 2 cos <t 2 , 7 2 sin <r 2 ) for the three orbits in Fig. 1. The plane 
is defined as <7i = 270° ,&i > 0. To obtain these, we adopt 
the initial conditions of orbits from Fig. 1, and then integrate 
the corresponding Hamiltonian equations based on the ex- 
pansion as Eq. (2,3). 

Fig. 9c indicates case c has a chaotic orbit. On the con- 
trary, the invariant curves for cases a and b (Fig. 9a,b) imply 
that both these two orbits are regular. That is, no matter 
whether the apsidal corotation happens or not, the system 
trapped in the 3:1 MMR could be stable. From this point of 
view we may argue the apsidal corotation only has a limited 
contribution to the stability of the system. 

Because the Hamiltonian is expressed explicitly in dj, 
and the integrals of motion J sum , Jdif affect the system as 
constrains, we can only define an energy surface (H = Ho) 
with certain pre-determined J SU m,./dif, then calculate the 
surface of section. We find the structure of the phase space 
and the stability of the system depends sensitively not only 
on H but also on J sum and Jdif. For example, with the same 
values of H, J aU m and Jdif as the ones in Fig. lb, we calculate 
tens of orbits with different initial conditions. On the surface 
of section, these orbits occupy a big regular-motion region. 
This result does not conflict with the fact that there are only 
a small fraction of stable systems in numerical simulations, 
in fact, when we change H, (or J sum , Jdif) just a little, we 
find remarkable chaotic region on the surface of section. The 
application of the surface-of-section technique is limited by 
the high dimension of the system. This will be discussed in 
detail in our future paper. 



5 SUMMARY AND DISCUSSION 

With hundreds of numerical simulations of the planetary 
system of the 55 Cancri, we find the third planet has a 
very weak influence on the motion of the inner two plan- 
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Figure 9. The projection on the surfaces of sections for case a, b and c in Fig. 1. 



ets. We confirm the inner two could be trapped in a 3:1 
mean motion resonance and three different types of motion 
are found. Judging from the Lyapunov character indicators 
and the surviving time of integration, two of them (case a, 
b) are practically stable, so that the real system could be 
running in one of these configurations. 

Via a new Hamiltonian expansion which is suitable for 
high-eccentricities planar three-body problem, we study the 
dynamics of the different configurations. We discuss the vari- 
ations of eccentricities and resonant angles, explain the hap- 
penings of different evolving types, and give a criterion of 
the occurrence of the apsidal corotation. 

The surfaces of section for the three types of motion are 
calculated and they reveal the stabilities of systems with or 
without the apsidal corotation. With these results we argue 
that the stability of the system is mainly due to the 3:1 
MMR, and the apsidal corotation has only a limited contri- 
bution. We would like to mention that this method can also 
be applied to other extra-solar planetary systems with other 
mean motion resonances. 

Numerical simulations suggest at least ~ 10 percent of 
systems are in a 3:1 MMR, while the Hamiltonian analy- 
ses give an upper limit of ~ 30%, which however will drop 
down after considering the stability. We also list the initial 
conditions leading to this MMR. 

Systems with different motion configurations have dif- 
ferent energy (H) levels. The apsidal corotation happens 
when the Hamiltonian approaches the extreme value. So, 
if the two planets are captured into current configuration 
through orbital migration caused by the action of non- 
conservative forces, the system should have an extremum of 
energy so that the apsidal corotation happens (Kley 2003). 
When the future observations would reveal more accurate 
properties of this system, we may consider what signatures 
of the migration are still presented in this system. 

The masses of planets adopted in this paper are the 
values from the orbital solutions when assuming sini = 1. 
As for the situations of sini < 1, our initial analysis with 
the Hamiltonian get the very similar results in a wide range 
of sini. This is also consistent with the results in (Beauge 
et al 2003). We have also analysed the possible motion con- 
figurations and their stabilities if the initial eccentricities 
differ from the values adopted above. With the help of the 



Hamiltonian model, we find that a higher e2 favors the es- 
tablishment of a 3:1 MMR. Anyway, we'd like to leave more 
details of these to our future paper. 

Last but not least, the general relativity effect may af- 
fects the secular dynamics of the 55 Cancri system, since the 
Companion B and C are quite close to the central star. For 
the innermost planet, the orbital precession caused by the 
general relativity effect, is calculated. Although this perias- 
tron shift, ~ 1.66 x 10 -6 radians per orbit, is quite small, 
it's about three times larger than that of Mercury in our 
Solar system. 
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